#ifndef _MESH_
#define _MESH_

#include "Domain.h"
#include "Dofs.h"
#include "Grid.h"
#include <map>
#include <functional>
#include <Eigen/Sparse>
#include <iostream>
#include <algorithm>
#include <string>
#include <fstream>
#include <sstream>
using namespace std;
/**
 * @brief 定义域划分成网格Mesh,Mesh包含小网格Grid,每个Grid的类型由TemplateElement确定
 * 
 * @tparam int DIM 维度
 */
template <unsigned int DIM>
class Mesh
{
public:
    /**
     * @brief 默认构造函数
     */
    Mesh() = default;
    
    /**
     * @brief 默认析构函数
     */
    ~Mesh() = default;
    
    /**
     * @brief 返回自由度的总个数
     * @details _GridList.size()*getNumDofs();
     * @return 
     */
    virtual int getTotalNumDofs(int _order_of_element) = 0;
    /**
     * @brief 返回网格的总个数
     * @details _GridList.size()
     * @return [description]
     */
    virtual int getNumGrid() = 0;
    /**
     * @brief 返回存储Grid的vector
     * @details _GridList
     * @return [description]
     */
    virtual vector<Grid<DIM>>& getGridVector() = 0;
    /**
     * @brief 返回序号为_index的Grid
     * 
     * @param _index Grid的序号
     * @return _GridList[i];
     */
    virtual Grid<DIM>& getGridVector(int _index) = 0;
    /**
     * @brief   返回vector，存储边界上的自由度
     * @details 对所有grid,依据seg找出边界的grid,再确定对应边界的自由度下标是哪几个
     *
     * @param component 试图对边界分段，并标号
     * @return 存储边界上自由度的vector
     */
    virtual vector<Dofs<DIM>>& getBoundaryDofs(int _order_of_element, int component = 0) = 0;

protected:

};

template <unsigned int DIM>
class StructedMesh : public Mesh<DIM>
{
public:
    StructedMesh() = default;
    /**
    * @brief 含参构造函数,根据给定的定义域和划分段数来划分网格。
    * @param domain 定义域
    * @param seg 每个维度要划分的段数
    */
    StructedMesh(Domain<DIM> domain, vector<int> _seg);
    ~StructedMesh() = default;
   
    int getTotalNumDofs(int _order_of_element);
    int getNumGrid();
    vector<Grid<DIM>>& getGridVector();
    Grid<DIM>& getGridVector(int _index);
    vector<Dofs<DIM>>& getBoundaryDofs(int _order_of_element, int component = 0);
    /**
     * @brief 得到分段
     * 
     */
    const vector<int>& getSegment() const;
    /**
     * @brief 得到第i个gird的所有自由度信息
     * 
     */
    const vector<Dofs<DIM>>& getDofsList(int i) const;
    /**
     * @brief 给定模板单元的阶数，在每个Grid上划定自由度。
     * @details 仅实现二维一阶、二阶模板单元
     * 
     * @param _order_of_element 模板单元的阶数
     */
    void distributeDofs(int _order_of_element);
protected:
    vector<int> seg;/*< 存放各维度分段数*/
    vector<vector<Dofs<DIM>>> DofsList;
    vector<Grid<DIM>> gridList;/*< 存放所有的Grid*/
    vector<Dofs<DIM>> boundaryDofs;/*< 存放所有边界上的自由度编号*/ 
};

template <unsigned int DIM>
StructedMesh<DIM>::StructedMesh(Domain<DIM> domain, vector<int> _seg)
{
    assert(DIM==2 && _seg.size() == DIM);
	Point<2> p0 = domain.getVertex(0);
	Point<2> p2 = domain.getVertex(2);
	double xLen = (p2[0] - p0[0])/(1.0*_seg[0]);
	double yLen = (p2[1] - p0[1])/(1.0*_seg[1]);
	for(int j = 0; j < _seg[1]; ++j)
	{
		for(int i = 0; i < _seg[0]; ++i)
		{
			vector<Point<2>> pnts;
			pnts.push_back({p0[0]+i*xLen,p0[1]+j*yLen});
			pnts.push_back({p0[0]+(i+1)*xLen,p0[1]+j*yLen});
			pnts.push_back({p0[0]+(i+1)*xLen,p0[1]+(j+1)*yLen});
			pnts.push_back({p0[0]+i*xLen,p0[1]+(j+1)*yLen});
			Grid<2> grid(pnts);
			gridList.push_back(grid);
		}
	}
	seg = _seg;
};

template <unsigned int DIM>
const vector<int>& StructedMesh<DIM>::getSegment() const
{
    return seg;
};

template <unsigned int DIM>
const vector<Dofs<DIM>>& StructedMesh<DIM>::getDofsList(int i) const
{
    return DofsList[i];
};

template <unsigned int DIM>
int StructedMesh<DIM>::getNumGrid()
{
    return gridList.size();
};
template <unsigned int DIM>
vector<Grid<DIM>>& StructedMesh<DIM>::getGridVector()
{
    return gridList;
};
template <unsigned int DIM>
Grid<DIM>& StructedMesh<DIM>::getGridVector(int _index)
{
    return gridList[_index];
};

template <unsigned int DIM>
vector<Dofs<DIM>>& StructedMesh<DIM>::getBoundaryDofs(int _order_of_element, int component)
{
    assert(DIM ==2);
    if(_order_of_element == 1)
    {
    	for(int k = 0; k < seg[0]; k++)
            boundaryDofs.push_back(DofsList[k][0]);
    	for(int k = 0; k <= seg[0] * (seg[1] - 1); k = k + seg[0])
            boundaryDofs.push_back(DofsList[k][3]);
    	for(int k = seg[0] - 1; k < seg[0] * seg[1]; k = k + seg[0])
            boundaryDofs.push_back(DofsList[k][1]);
    	for(int k = seg[0] * (seg[1] - 1); k < seg[0] * seg[1]; k++)
            boundaryDofs.push_back(DofsList[k][2]);
    }
    else if(_order_of_element == 2)
    {
    	for(int k = 0; k < seg[0]; k++){
            boundaryDofs.push_back(DofsList[k][0]);
            boundaryDofs.push_back(DofsList[k][1]);}
    	for(int k = 0; k <= seg[0] * (seg[1] - 1); k = k + seg[0]){
            boundaryDofs.push_back(DofsList[k][7]);
            boundaryDofs.push_back(DofsList[k][6]);}
    	for(int k = seg[0] - 1; k < seg[0] * seg[1]; k = k + seg[0]){
            boundaryDofs.push_back(DofsList[k][2]);
            boundaryDofs.push_back(DofsList[k][3]);}
    	for(int k = seg[0] * (seg[1] - 1); k < seg[0] * seg[1]; k++){
            boundaryDofs.push_back(DofsList[k][5]);
            boundaryDofs.push_back(DofsList[k][4]);}
    }
    else
    	abort(3);
    return boundaryDofs;
};


//假定所有单元类型、阶数都相同
template <unsigned int DIM>
void StructedMesh<DIM>::distributeDofs(int _order_of_element)
{
    //以下示范DIM=2，_order_of_element = 1 | 2;
    assert(DIM ==2);
    int k = 0;
    vector<Dofs<2>> _dofs;
    if(_order_of_element == 1)
    {
        for (Grid<DIM> &grid : gridList)
        {
            _dofs.clear();
            vector<Point<DIM>> p = grid.getVertices();
            _dofs.push_back({{p[0][0], p[0][1]}, k + k / seg[0]});
            _dofs.push_back({{p[1][0], p[1][1]}, k + k / seg[0] + 1});
            _dofs.push_back({{p[2][0], p[2][1]}, k + k / seg[0] + seg[0] + 2});
            _dofs.push_back({{p[3][0], p[3][1]}, k + k / seg[0] + seg[0] + 1});
            
            DofsList.push_back(_dofs);
            k++;
        }

    }
    else if(_order_of_element == 2)
    {
        vector<Dofs<DIM>> _dofs;
        for (Grid<DIM> &grid : gridList)
        {
            _dofs.clear();
            int a = seg[0];
            vector<Point<DIM>> p = grid.getVertices();
            _dofs.push_back({{p[0][0], p[0][1]}, k / a * (4 * a + 2) + k % a * 2});
            _dofs.push_back({{(p[0][0] + p[1][0])/2.0, (p[0][1] + p[1][1])/2.0}, k / a * (4 * a + 2) + k % a * 2 + 1});
            _dofs.push_back({{p[1][0], p[1][1]}, k / a * (4 * a + 2) + k % a * 2 + 2});
            _dofs.push_back({{(p[1][0]+p[2][0])/2.0, (p[1][1] + p[2][1])/2.0}, k / a * (4 * a + 2) + 2 * a + 1 + k % a * 2 + 2});
            _dofs.push_back({{p[2][0],p[2][1]}, (k / a + 1) * (4 * a + 2) + k % a * 2 + 2});
            _dofs.push_back({{(p[2][0]+p[3][0])/2.0, (p[2][1] + p[3][1])/2.0}, (k / a + 1) * (4 * a + 2) + k % a * 2 + 1});
            _dofs.push_back({{p[3][0],p[3][1]}, (k / a + 1) * (4 * a + 2) + k % a * 2});
            _dofs.push_back({{(p[3][0]+p[0][0])/2.0, (p[3][1]+p[0][1])/2.0}, k / a * (4 * a + 2) + 2 * a + 1 + k % a * 2});
            _dofs.push_back({{(p[0][0]+p[2][0])/2.0, (p[0][1]+p[2][1])/2.0}, k / a * (4 * a + 2) + 2 * a + 1 + k % a * 2 + 1});

            DofsList.push_back(_dofs);
            k++;
        }
    }
    else
        abort(3); 
};

template <unsigned int DIM>
int StructedMesh<DIM>::getTotalNumDofs(int _order_of_element)
{
    if(_order_of_element == 1)
    	return (seg[0] + 1) * (seg[1] + 1);
    else if(_order_of_element == 2)
    	return (2 * seg[0] + 1) * (2 * seg[1] + 1);
    else
    	abort(3);
};


template<unsigned int DIM>
class UnstructedMesh: public Mesh<DIM>
{
public:

	virtual void InputDofData() = 0;
    	virtual void InputEdgeData() = 0;
   	virtual void InputEleData() = 0;

};

class EasyMesh : public  UnstructedMesh<2>
{
private:
    string _MeshName;
    long int _ndofs;
    long int _nele;
    long int _nbndedge;
    vector<Dofs<2>> _DofsList;
    vector<Dofs<2>> _BndDofs;
    vector<vector<int>> _ElementData;
    vector<vector<int>> _BndEdges; 

public:
    EasyMesh();
    EasyMesh(string meshname);
    void InputDofData();
    void InputEdgeData();
    void InputEleData();
    int getNumGrid(){return _nele;};
    int getTotalNumDofs(){return _ndofs;};
    int getBndedgeNum(){return _nbndedge;};
    Dofs<2> getDofs(int i){return _DofsList[i];};
    vector<int> getElement(int i){return _ElementData[i];};
    vector<Dofs<2>>& getBoundary(){return _BndDofs;};
    //bool IsBoundary(int index){return !(std::find(_BndDofs.begin(),_BndDofs.end(),index) == _BndDofs.end());};
    vector<int> getBndEdge(int i){return _BndEdges[i];};    
};

EasyMesh::EasyMesh(string meshname)
{
    _MeshName = meshname;
    InputDofData();
    InputEdgeData();
    InputEleData();
}

void EasyMesh::InputDofData()
{
    ifstream fin;
    fin.open(_MeshName+".n");
    string line;
    stringstream ss;
    getline(fin, line);
    ss.clear();
    ss.str(line);
    ss >> _ndofs;
    _DofsList.resize(_ndofs);

    int a, bndflag;
    double x, y;
    for (int i = 0; i < _ndofs; i++)
    {
	getline(fin, line);
	ss.clear();
	ss.str(line);
	ss >> a >> x >> y >> bndflag;

	Point<2> temp_point({x, y});
	Dofs<2> temp_dof(temp_point, i);
	_DofsList[i] = temp_dof;
	if (bndflag)
	    _BndDofs.push_back(temp_dof);	    
    }
}

void EasyMesh::InputEdgeData()
{
    ifstream fin;
    fin.open(_MeshName+".s");
    string line;
    stringstream ss;
    getline(fin, line);
    ss.clear();
    ss.str(line);
    int _nedge;
    ss >> _nedge;
    int a, d0, d1, bndflag;
    for (int i = 0; i < _nedge; i++)
    {
	getline(fin, line);
	ss.clear();
	ss.str(line);
	ss >> a >> d0 >> d1 >> a >> a >> bndflag;
	if (bndflag)
	{
	    vector<int> temp_edge(2);
	    temp_edge[0] = d0;
	    temp_edge[1] = d1;
	    _BndEdges.push_back(temp_edge);
 	}
    }
    _nbndedge = _BndEdges.size();
}

void EasyMesh::InputEleData()
{
    ifstream fin;
    fin.open(_MeshName+".e");
    string line;
    stringstream ss;
    getline(fin, line);
    ss.clear();
    ss.str(line);
    ss >> _nele;
    _ElementData.resize(_nele);

    int a, d0, d1 ,d2;
    for (int i = 0; i < _nele; i++)
    {
	getline(fin, line);
	ss.clear();
	ss.str(line);
	ss >> a >> d0 >> d1 >> d2;
	_ElementData[i].resize(3);
	_ElementData[i][0] = d0;
	_ElementData[i][1] = d1;
	_ElementData[i][2] = d2;
    }
}

#else
#endif
